Introduction

k-Fold Cross Validation

Data:

For this CS example, we use the common dataset bwght from Wooldridge: Introductory econometrics.

##  [1] "bwght"    "cigtax"   "cigprice" "faminc"   "fatheduc" "motheduc"
##  [7] "parity"   "male"     "white"    "cigs"     "bwghtlbs" "packs"

Basic data plot:

plot(BWGHT[,c(1,4:7)])


Three simple linear models to be used for comparison:

# Model estimation: three alternative specifications
Model.1 <- lm(bwght ~ cigs+faminc, data=BWGHT)
Model.2 <- lm(bwght ~ cigs+I(cigs^2)+faminc, data=BWGHT)
Model.3 <- lm(bwght ~ cigs+faminc+parity+male, data=BWGHT)

Regression table outputs

summary(Model.1)
## 
## Call:
## lm(formula = bwght ~ cigs + faminc, data = BWGHT)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -96.061 -11.543   0.638  13.126 150.083 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 116.97413    1.04898 111.512  < 2e-16 ***
## cigs         -0.46341    0.09158  -5.060 4.75e-07 ***
## faminc        0.09276    0.02919   3.178  0.00151 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 20.06 on 1385 degrees of freedom
## Multiple R-squared:  0.0298, Adjusted R-squared:  0.0284 
## F-statistic: 21.27 on 2 and 1385 DF,  p-value: 7.942e-10
summary(Model.2)
## 
## Call:
## lm(formula = bwght ~ cigs + I(cigs^2) + faminc, data = BWGHT)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -96.243 -11.789   0.895  13.131 149.958 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 117.218964   1.058478 110.743  < 2e-16 ***
## cigs         -0.776637   0.208464  -3.726 0.000203 ***
## I(cigs^2)     0.012256   0.007329   1.672 0.094686 .  
## faminc        0.089961   0.029217   3.079 0.002117 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 20.05 on 1384 degrees of freedom
## Multiple R-squared:  0.03176,    Adjusted R-squared:  0.02966 
## F-statistic: 15.13 on 3 and 1384 DF,  p-value: 1.078e-09
summary(Model.3)
## 
## Call:
## lm(formula = bwght ~ cigs + faminc + parity + male, data = BWGHT)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -96.498 -11.834   0.854  13.322 152.621 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 112.38966    1.59085  70.647  < 2e-16 ***
## cigs         -0.47501    0.09127  -5.205 2.24e-07 ***
## faminc        0.10219    0.02914   3.507 0.000468 ***
## parity        1.64607    0.60237   2.733 0.006363 ** 
## male          3.16310    1.07404   2.945 0.003283 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.96 on 1383 degrees of freedom
## Multiple R-squared:  0.04082,    Adjusted R-squared:  0.03804 
## F-statistic: 14.71 on 4 and 1383 DF,  p-value: 8.935e-12

kFCV

Use ?cvFit from the {cvTools} to familiarize with arguments to the cvFit function. The following settings are used here:

  • 5-folds for cross-validation,
  • 100 repetitions of the whole kFCV procedure (i.e. variance and distribution of the \(CV(k)\) statistic can be empirically evaluated),
  • default evaluation: using root mean squared prediction error: rtmspe
Model 1
cv.Model.1 <- cvFit(Model.1, data = BWGHT, y = BWGHT$bwght, cost = rtmspe, 
                    K = 5, R = 100, seed = 1234)
cv.Model.1
## 5-fold CV results:
##       CV 
## 11.10385
Model 2
cv.Model.2 <- cvFit(Model.2, data = BWGHT, y = BWGHT$bwght, cost = rtmspe, 
                    K = 5, R = 100, seed = 1234)
cv.Model.2
## 5-fold CV results:
##       CV 
## 11.11938
Model 3
cv.Model.3 <- cvFit(Model.3, data = BWGHT, y = BWGHT$bwght, cost = rtmspe, 
                    K = 5, R = 100, seed = 1234)
cv.Model.3
## 5-fold CV results:
##       CV 
## 11.05869

Evaluation table and evaluation plot for k-fold CV output

cvFits <- cvSelect(Model_1 = cv.Model.1, Model_2 = cv.Model.2, 
                   Model_3 = cv.Model.3)
cvFits
## 
## 5-fold CV results:
##       Fit       CV
## 1 Model_1 11.10385
## 2 Model_2 11.11938
## 3 Model_3 11.05869
## 
## Best model:
##        CV 
## "Model_3"
summary(cvFits)
## 
## 5-fold CV results:
##       Fit  CV.Min. CV.1st Qu. CV.Median  CV.Mean CV.3rd Qu.  CV.Max.
## 1 Model_1 11.07549   11.09232  11.10153 11.10385   11.11035 11.16312
## 2 Model_2 11.08358   11.10542  11.11707 11.11938   11.13204 11.17687
## 3 Model_3 11.02297   11.04122  11.05780 11.05869   11.07158 11.11977
## 
## Best model:
##        CV 
## "Model_3"
# different plot outputs are available - uncomment to display
#   plot(cvFits, method = "bw")
#   plot(cvFits, method = "xy")
#   plot(cvFits, method = "dot")
plot(cvFits, method = "density") 

Additional topics